//************************* BEZIER CURVES *****************************************
pt cubicBezier(pt A, pt B, pt C, pt D, float t) {
  return( s( s( s(A, t, B), t, s(B, t, C) ), t, s( s(B, t, C), t, s(C, t, D) ) ) );
}
void drawCubicBezier(pt A, pt B, pt C, pt D) { 
  beginShape();  
  for (float t=0; t<=1; t+=0.02) {
    cubicBezier(A, B, C, D, t).v();
  };  
  endShape();
}

void showArcThrough (pt A, pt B, pt C) {
  if (abs(dot(V(A, B), R(V(A, C))))<0.01*d2(A, C)) {
    show(A, C); 
    return;
  }
  pt O = CircumCenter ( A, B, C); 
  float r=d(O, A);
  vec OA=V(O, A), OB=V(O, B), OC=V(O, C);
  float b = angle(OA, OB), c = angle(OA, OC); 
  if (0<c && c<b || b<0 && 0<c)  c-=TWO_PI; 
  else if (b<c && c<0 || c<0 && 0<b)  c+=TWO_PI; 
  beginShape(); 
  v(A); 
  for (float t=0; t<1; t+=0.01) v(R(A, t*c, O)); 
  v(C); 
  endShape();
}
//************************* CIRCLES *****************************************
float interpolateAngles(float a, float t, float b) {
  if (b<a) b+=TWO_PI; 
  float m=(t-a)/(b-a); 
  if (m>PI) m-=TWO_PI; 
  return m;
}


pt ptOnCircle(float a, pt O, float r) {
  return P(r*cos(a)+O.x, r*sin(a)+O.y);
} 
//************************* EDGES *****************************************
boolean edgesIntersect(pt A, pt B, pt C, pt D) {
  return cw(A, B, C)!=cw(A, B, D) && cw(C, D, A)!=cw(C, D, B) ;
}
boolean edgesIntersect(pt A, pt B, pt C, pt D, float e) {
  return ((A.isLeftOf(C, D, e) && B.isLeftOf(D, C, e))||(B.isLeftOf(C, D, e) && A.isLeftOf(D, C, e)))&&
    ((C.isLeftOf(A, B, e) && D.isLeftOf(B, A, e))||(D.isLeftOf(A, B, e) && C.isLeftOf(B, A, e))) ;
}
pt linesIntersection(pt A, pt B, pt C, pt D) {
  vec AB = V(A, B);  
  vec CD = V(C, D);  
  vec N=R(CD); 
  vec AC = V(A, C); 
  float s = dot(AC, N)/dot(AB, N); 
  return T(A, s, AB);
}

//************************* interpolation  *****************************************
float quadraticInterpolant(float t, float A, float B, float C) {
  return (t-1)*(t-2)*A/2 + t*(2.-t)*B + t*(t-1)*C/2;
}
pt PtOnParabola(float t, pt A, pt B, pt C) {
  return P( quadraticInterpolant(2*t, A.x, B.x, C.x), quadraticInterpolant(2*t, A.y, B.y, C.y) );
}
pt MidCoursePtOnCubic(pt A, pt B, pt C, pt D) {
  return A(PtOnParabola(3./4, A, B, C), PtOnParabola(1./4, B, C, D));
}
pt B(pt A, pt B, pt C, float s) {
  return( s(s(B, s/4., A), 0.5, s(B, s/4., C)));
};                          // returns a tucked B towards its neighbors
pt F(pt A, pt B, pt C, pt D, float s) {
  return( s(s(A, 1.+(1.-s)/8., B), 0.5, s(D, 1.+(1.-s)/8., C)));
};    // returns a bulged mid-edge point 
pt s(pt A, float s, pt B) {
  return(new pt(A.x+s*(B.x-A.x), A.y+s*(B.y-A.y)));
};


//************************* spirals  *****************************************
void drawSpiral(pt A, pt B, float t, pt C, pt D) { 
  float a = spiralAngle(A, B, C, D); 
  float s = spiralScale(A, B, C, D);
  pt G = spiralCenter(a, s, A, C); 
  show(G);
  pt P = spiralPt(A, G, s, a, t); 
  pt Q = spiralPt(B, G, s, a, t);
  show(P, Q);
}
pt spiralPt(pt A, pt G, float s, float a) {
  return L(G, s, R(A, a, G));
}  
pt spiralPt(pt A, pt G, float s, float a, float t) {
  return L(G, pow(s, t), R(A, t*a, G));
} 
//pt spiralPt(pt A, pt G, float s, float a, float t) {return L(G,(1.-t)+t*s,R(A,t*a,G));} 
pt spiralCenter(pt A, pt B, pt C, pt D) { // computes center of spiral that takes A to C and B to D
  float a = spiralAngle(A, B, C, D); 
  float z = spiralScale(A, B, C, D);
  return spiralCenter(a, z, A, C);
}
float spiralAngle(pt A, pt B, pt C, pt D) {
  return angle(V(A, B), V(C, D));
}
float spiralScale(pt A, pt B, pt C, pt D) {
  return d(C, D)/d(A, B);
}
pt spiralCenter(float a, float z, pt A, pt C) {
  float c=cos(a), s=sin(a);
  float D = sq(c*z-1)+sq(s*z);
  float ex = c*z*A.x - C.x - s*z*A.y;
  float ey = c*z*A.y - C.y + s*z*A.x;
  float x=(ex*(c*z-1) + ey*s*z) / D;
  float y=(ey*(c*z-1) - ex*s*z) / D;
  return P(x, y);
}
/************************************* triangles **************************************/
// measure
boolean cw(pt B, pt C) {
  return C.y*B.x>C.x*B.y;
}  // cw(B,C): Boolean true if smallest angle turn from OB to OC is cloclwise (cw), where O=(0,0)
boolean cw(pt A, pt B, pt C) {
  return (C.y-A.y)*(B.x-A.x)>(C.x-A.x)*(B.y-A.y);
} // cw(A,B,C): Boolean true if A-B-C makes a clockwise turn at B
float area(pt A, pt B, pt C) {
  return 0.5*((C.y-A.y)*(B.x-A.x)+(A.x-C.x)*(B.y-A.y));
} // area(A,B,C): signed Triangle area (positive if A-B-C is clockwise
float x(pt A, pt B, pt C) {
  return ((C.x-A.x)*(B.x-A.x)+(C.y-A.y)*(B.y-A.y))/(sq(C.x-A.x)+sq(C.y-A.y));
} // x(A,B,C): local x-coordinate of B in system (AC,R(AC),A)
float y(pt A, pt B, pt C) {
  return ((C.y-A.y)*(B.x-A.x)+(A.x-C.x)*(B.y-A.y))/(sq(C.x-A.x)+sq(C.y-A.y));
} // y(A,B,C): local y-coordinate of B in system (AC,R(AC),A)
float ay(pt A, pt B, pt C) {
  return ((C.y-A.y)*(B.x-A.x)+(A.x-C.x)*(B.y-A.y))/sqrt(sq(C.x-A.x)+sq(C.y-A.y));
} // ay(A,B,C): signed distance from A to the projection of B onto line(A,R(AC))
pt projection(pt A, pt B, pt C) {
  return L(A, y(A, B, C), C);
} // returns projection of B on line(A,C)
float thickness(pt A, pt B, pt C) {
  return min(abs(ay(A, B, C)), abs(ay(B, C, A)), abs(ay(C, A, B)));
} 

// area , moment
float triangleArea(pt A, pt B, pt C) {
  return(dot(R(V(A, B)), V(A, C))/2.);
};
float triangleMoment(pt A, pt B, pt C) {
  float b = A.disTo(B); 
  vec T=A.makeVecTo(B); 
  T.normalize();
  vec N = R(T); 
  vec AC=V(A, C); 
  float h = dot(AC, N);
  float a = dot(AC, T);
  return ( b*b*b*h - a*b*b*h + a*a*b*h + b*h*h*h )/36.;
};

// render
void show(pt A, pt B, pt C, float r) {
  if (thickness(A, B, C)<abs(2*r)) return;
  float s=r; 
  if (cw(A, B, C)) s=-s; 
  beginShape();  
  Offset(A, B, C, s).v(); 
  Offset(B, C, A, s).v(); 
  Offset(C, A, B, s).v(); 
  endShape(CLOSE);
}
void showAltitude(pt A, pt B, pt C) {
  vec V=R(V(A, C)); 
  vec U=S(y(A, B, C), V); 
  arrow(B, U);
}

// centers and their radii: http://www.jimloy.com/geometry/centers.htm
pt CircumCenter (pt A, pt B, pt C) {
  vec AB = V(A, B); 
  vec AC = R(V(A, C)); 
  return T(A, 1./2/dot(AB, AC), W(-n2(AC), R(AB), n2(AB), AC));
}; // CircumCenter(A,B,C): center of circumscribing circle, where medians meet)
float circumRadius (pt A, pt B, pt C) {
  float a=d(B, C), b=d(C, A), c=d(A, B), s=(a+b+c)/2, d=sqrt(s*(s-a)*(s-b)*(s-c)); 
  return a*b*c/4/d;
} // radiusCircum(A,B,C): radius of circumcenter

pt OrthoCenter(pt A, pt B, pt C) {
  float a=d2(B, C), b=d2(C, A), c=d2(A, B), x=(a+b-c)*(a-b+c), y=(a+b-c)*(-a+b+c), z=(a-b+c)*(-a+b+c), t=x+y+z; 
  x=x/t; 
  y=y/t; 
  z=z/t; 
  return W(x, A, y, B, z, C);
} // OrthoCenter(A,B,C): intersection of altitudes
pt SpiekerCenter (pt A, pt B, pt C) {
  float ab=d(A, B), bc=d(B, C), ca=d(C, A), s=ab+bc+ca; 
  return W(ab/s, A(A, B), bc/s, A(B, C), ca/s, A(C, A));
}    // SpiekerCenter(A,B,C): center of mass of the perimeter 
float spiekerRadius (pt A, pt B, pt C) {
  float ab=d(A, B), bc=d(B, C), ca=d(C, A), s=ab+bc+ca; 
  return s/(2*PI);
}    // spiekerRadius(A,B,C): radius of circle with same perimeter
pt InCenter (pt A, pt B, pt C) {
  float Z=area(A, B, C), a=B.disTo(C), b=C.disTo(A), c=A.disTo(B), s=a+b+c, r=2*Z/s, R=a*b*c/(2*r*s); 
  return W(a/s, A, b/s, B, c/s, C);
} // InCenter(A,B,C): incenter (center of inscribed circle)
float inRadius (pt A, pt B, pt C) {
  float Z=area(A, B, C), a=d(B, C), b=d(C, A), c=d(A, B), s=a+b+c;  
  return 2*Z/s;
} //inRadius(A,B,C): radius of incircle  

float radiusMonotonic (pt A, pt B, pt C) {                                                  // size of bubble pushed through (A,C) and touching B, >0 when ABC is clockwise
  float a=d(B, C), b=d(C, A), c=d(A, B);
  float s=(a+b+c)/2; 
  float d=sqrt(s*(s-a)*(s-b)*(s-c)); 
  float r=a*b*c/4/d;
  if (abs(angle(A, B, C))>PI/2) r=sq(d(A, C)/2)/r;
  if (abs(angle(C, A, B))>PI/2) r=sq(d(C, B)/2)/r;
  if (abs(angle(B, C, A))>PI/2) r=sq(d(B, A)/2)/r;
  if (cw(A, B, C)) r=-r;
  return r;
};

pt Offset(pt A, pt B, pt C, float r) {  // intersection of r-offsets of lines (A,B) and (B,C)
  float a = angle(V(B, A), V(B, C))/2;
  float d = r/sin(a); 
  vec N = U(A(U(V(B, A)), U(V(B, C)))); 
  return T(B, d, N);
};

pt Shadow(pt A, pt B, pt C) {
  return L(A, x(A, B, C), C);
};  // Projection of B on line (B,C)





